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COMPUTER PROGRAM TO CALCULATE CONVECTIVE HEAT -TRANSFER 
COEFFICIENTS FROM HARMONIC TEMPERATURE OSCILLATIONS 
by Geraldine E. Amling and Ronald G. Huff 
Lewis Research Center 

SUMMARY 

Two FORTRAN IV version 13 programs are described which obtain the analytical 
solutions for convective heat-transfer coefficients related to either the phase lag or the 
amplitude ratio. The first, entitled BL, uses the ratio of the wall to driven hot fluid 
temperature to calculate the coefficients. The second, entitled JCL, uses the phase-lag 
angle between the sinusoidally driven fluid and wall temperature to calculate the heating 
fluid and coolant convective coefficients. Criteria for the choice of either of the two 
methods and their mathematical formulation are described in detail in reference 1. 

INTRODUCTION 

The analytical expressions to determine the convective heat-transfer coefficients for 
a wall separating two fluids, one having an oscillating temperature, were derived in ref- 
erence 1 from a second-order partial differential equation. These equations lead to 
computational procedures that are iterative and, therefore, time consuming and tedious. 
To alleviate this difficulty, two computer programs, BL and JCL, were written to pro- 
vide the means for obtaining the numerical solutions to either the amplitude ratio or the 
phase-lag angle method quickly and accurately. 

The data published in reference 1 were obtained by the use of these programs, and 
their descriptions and listings, which are shown herein, should prove helpful to anyone 
desiring to investigate further. 

Specifically, the BL program calculates the convective heat-transfer coefficients 
related to the amplitude ratio for a plate separating two fluids, one of the fluids having 
an oscillating temperature. The JCL program, on the other hand, calculates the heat- 
transfer coefficients related to the phase lag. 



STATEMENT OF THE PROBLEM 


Figure 1(a) shows a wall that separates two moving fluids. The hot-gas-side fluid 
temperature Tq is greater than the coolant temperature T c , causing heat to flow 
through the wall in the positive x direction. The problem is to determine the convec- 
tive heat-transfer coefficients if either the hot-gas or coolant temperature is varied 
sinusoidally and if the temperature response of the wall is measured at only one point. 
Figure 1(b) shows what the wall temperature" might look like at any given instant in time 



(a) Wall temperature distribution shown at time zero before start of hot -gas temperature 
oscillation. 
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(b) Transient temperature response of plate to sinusoidally forced fluid temperature. 
Figure 1. - Basic heat-transfer model. 



METHOD OF ATTACK 


If the hot-gas temperature is made to vary sinusoidally, the wall temperature will 
respond sinusoidally but will lag the driven temperature by an angle cp. In addition, the 
amplitude of the wall temperature oscillation will be less than that of the driven tempera- 
ture. Both the phase lag cp and the ratio of the amplitude of the wall temperature to the 
amplitude of the driven temperature 0/AT are, among other things, functions of the 
convective heat-transfer coefficients h. Finding the relation between the convective 
heat-transfer coefficients and the phase lag cp, or the amplitude ratio 0^/AT, is re- 
quired in order to solve for the coefficients. 

To find this relation, the temperature response of the wall to a sinusoidally driven 
fluid temperature is derived. From this solution and steady-state heat-transfer condi- 
tions (which relate the ratio of the convective heat-transfer coefficients to the ratio of the 
temperature drops between fluids and wall), the absolute values of the coefficients can be 
determined as functions of either phase lag or amplitude ratio. The quantities that must 
be known as a function of time are the hot-gas temperature T G , the wall temperature at 
any point x (T(x, r)) and the coolant temperature T c . From these quantities either the 
phase lag or the amplitude ratio can be determined as well as the ratio of the convective 
heat-transfer coefficients. With the ratio of the coefficients, the frequency of the tem- 
perature oscillation, the wall material properties, and either cp or 0 m /ATQ, the con- 
vective heat-transfer coefficients h G and h £ can be calculated. 

It is also possible to calculate the coefficients by oscillating the coolant temperature 
sinusoidally. This approach requires the measurement of the same quantities as when 
the hot fluid temperature is oscillated. 

The controlling differential equations and boundary conditions are given next. The 
one- dimensional transient heat- conduction equation 1 


S 2 T 

ax 2 


1 9T 

a 3 r 


( 1 ) 


is solved for the temperature distribution in a wall that has convective heat transfer over 
its two surfaces (fig. 1). (Symbols are defined in the appendix. ) The hot-gas side fluid 
temperature is driven sinusoidally and is given in equation form by 

T G = AT G e" iCl,T + T G (2) 


1 


All equations are numbered to correspond with those in reference 1. 
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where Tq is given. It is not necessary to know the numerical value of AT^ if the 
phase lag is used for calculating the coefficients because the phase determination is inde- 
pendent of the absolute value of the temperature oscillation. If the amplitude ratio is 
used for calculating the coefficients, the ATq must be computed from the change in 
chamber pressure. The boundary condition at the hot-gas side of the wall surface 
x/L = 0 equates the heat transferred by convection to that conducted away from the sur- 
face into the wall. This is given in equation form as 

h G [ T G - T(0, t)J = -K - ST ^ Q > T ) (3) 

The boundary condition at the coolant surface (x/L = 1) equates the heat conducted to the 
surface with the heat transferred convectively to the coolant. In equation form this is 
given by 


h c [T(L, t) - T c ] - (4) 

The solution for the wall temperature as a function of time and location within the 
wall consists of the sum of the transient 0(x, r) and steady- state T_(x) solutions. The 
steady-state solution is found by setting 3T/9 r equal to zero in equation (1) and assum- 
ing the wall temperature at x/L = 0 and 1 to be given as Tq W and T cw , respectively. 
This yields the usual linear temperature distribution with distance and can be written as 

< 5 > 

The transient solution 0(x, r) is found by assuming the usual product solution: one 
a function of time only F(t) and the other a function of the location within the wall meas- 
ured from the heated side of the wall X(x). Equations (1), (3), and (4) are modified to 
account for the change in variable from T(x, t) to 0(x, r). The details of the solution of 
equation (1) using its pertinent boundary conditions (eqs. (2) to (4)) are shown in appen- 
dix B of reference 1. 

The relation, derived from the steady- state (or mean) condition that requires the 
heat transferred convectively to the wall to equal the heat transferred convectively to the 
coolant, is expressed mathematically as 
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( 6 ) 


h c _ T g - t (°) 
h G T(L) - T c 


From the definition of 


h = 


* 


(V 


The ratio of convective heat-transfer coefficients in equation (6) is defined as R. Calcu- 
lating the ratio R using equation (7) gives 



(8) 


where R is also given by 

T„ - T(0) 

R — 

T(L) - T c 


(9) 


The following assumptions have been used to obtain the solution: 

(1) The heat flows through the plate in the x direction only (one -dimensional heat 
conduction). 

(2) The wall or plate properties (density, specific heat, and thermal conductivity) 
are constant. 

(3) The convective heat-transfer coefficients on both sides of the wall are constant. 

(4) The coolant temperature is constant. 

The details of the analytical solutions are given in reference 1. However, the fol- 
lowing equations are restated herein because they are used by the computer program to 
calculate the heat-transfer coefficients. 

The equation involving amplitude ratio is 


"m _ 1 


(ca - db)sin r|L — - (da + cb)cos P" ( X/,L 3 + Ufa - eb)eos t)L — - (ea + fb)sin 7)L —1 e ’1 L ( X A‘) 

L L I L LJ 


a 


.{[„ 


+< (ea + fb)cos 77 L — + (fa 
L 


- eb)sin tjL ije^x/L) - J^(ca -db)cos tjL i + (da + cb)sin „L £j e ^[2-(xA)]j‘ 


1/2 


(B37) 
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where 


a = 1 - (i^ G + »^ c ) + e 2?7L ^ G + + 2^ G ^ c )sin 277 L - (1 + + ^ G )cos 2 ??lJ (B28) 

b = (i^ G + - 2x Pq^ c ) + e 2 ^ L b G + + 2 ^q^ cos 2tjL + ( 1 + + ^G^ sin 2??L 1 



(B29) 

c = (1 + ty c )cos 2rjL - ip c sin 2rjL 

(B30) 

d = i// cos 2rjL + (1 + i/O sin 217 L 

(B31) 

e “ 1 “ ^c 

(B32) 

II 

O 

(B33) 


The equation involving phase lag angle is 


3 [C+R(C+D) 2 2A + R(A + B) E 

\]yl + l]S + lj/ + - 0 

c R(F + G) c 2R(F + G) c 2R(F + G) 


(B51) 


where 

A s -sin (cp + rjx) + e 2 ^^ 1 cos (cp + rjx - 2rjL ) + e 2 ^^ 2 ^ J-x ^sin((p - rjx) 

- e 2 ^^ x ^cos(<p - Tjx + 2t]L) (B39) 

B s cos (cp + rjx) - e 2r ^ sin(<p + rjx - 2 t]L) - e 2 ^ 2 ^ x ^cos (cp - rjx ) 

+ e 2? ^ x ^sin(<p - rjx + 2rjlt) (B40) 
C = sin(cp + rjx) + e 2 ^ L sin(<p + rjx - 2?jL) + e 2r, ( 2 ^~ , ~ x '* sin(cp - rjx) 

+ e 2? ^ x ^sin(<p - rjx + 2?7L) (B41) 

D s -cos (cp + tjx) + e 2 ^ L cos (cp + rjx - 2?jL) - e 2?7 ^ 2 '* J-x ^cos((p - ? 7 x) 

+ e 271 ^^ cos{cp - rjx + 2rih) (B42) 
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E s sin(<j? + tjx) - e 2? ^ sin(cp + tjx - 2tjL) + e^ 7 ^ 2 ^ x ) s j n (^ _ ^ x ) 

- e 2T7(L_x) sin(^ - tjx + 2tjL) (B43) 
F s cos (<p + ?]x) + e 2? ^ sin(cp + tjx - 2tjL) - e 2 ^^ 2 ^ 1 x ^cos(^? - tjx) 

- e 27? ( L_X W(<p - tjx + 2tjL) (B44) 
G = -sin (cp + tjx) - e 2 ^^ 1 eos(<j? + tjx - 2tjL) + e 2 ^^ 2 ^ x ^sin(<p - tjx) 

+ e 2T ^^ _x ^cos(^ - tjx + 2tjL) (B45) 


The solution of this cubic equation can be obtained using the classical approach. The 
positive, nonzero, root will yield the desired convective heat-transfer parameter i^ c . 
From equation (8), then 

= R^ c (B52) 

Also, from equations (6) and (9), the following relation must be satisfied: 


h c T G - T <°> 
h Q " T(L) - T c 


(B49) 


DESCRIPTION OF BLAND JCL 

The computer program, BL, was written to calculate the heat-transfer coefficients 
h c and h^ by use of values of the thermal conductivity K, the density of the wall mate- 
rial p, the specific heat of the wall material c, thickness of wall material L, the ratio 
of distance measured from the heated surface wall to the thickness of the wall x/L, the 
mean temperature of the hot gas T n , the wall temperature as a function of time at x = 0 
T(0, t), the mean temperature of the coolant T c , the frequency of temperature oscillation 
f, and the amplitude ratio 0 /ATq . 

The coolant wall temperature T cw is initially set to a positive perturbation value in 
subroutine BALAN for the coolant temperature T c . The ratio of the convective heat- 
transfer coefficients R, is then calculated, and with a linear relation of conductive heat 
transfer across the plate both heat-transfer coefficients h £ and h Q can be obtained. 
With these quantities, a calculated amplitude ratio is found and compared with the meas- 
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ured ratio. Then, by means of linear interpolation in subroutine BALAN, T cw is al- 
tered until the difference between the two ratios is less than 10 - ^. 

It should be noted that the 6 / ATq ratio as a function of the coolant wall temperature 
has a pole at T(x, t) (fig. 2), which imposes a heavy restriction on the allowable error 
for the value of the amplitude ratio. If the error in the amplitude ratio exceeds 10 per- 
cent of the true value, the iteration will converge to a value for the coolant wall temper- 
ature in excess of the gas wall temperature. 

The computer program JCL was written to calculate heat-transfer coefficients h 
and h^, by use of values of the thermal conductivity K, the density of wall material p, 
the specific heat of wall material c, thickness of wall material L, the ratio of distance 
measured from the heated surface of the wall to the thickness of the wall x/L, the mean 
temperature of the hot gas Tq, the wall temperature as a function of distance (measured 
into wall from hot surface) T(x, r), the mean temperature of the coolant T c , the frequen- 
cy of temperature oscillation f, and the phase lag angle 9 . The coolant wall tempera- 
ture T... is initially set to equal the measured temperature of the wall T(x, t). Thera- 
tio R, of the convective heat-transfer coefficients is calculated, and the coolant wall 
temperature value is corrected for the wall temperature drop due to conduction. Only the 
positive, nonzero root is used when solving the cubic equation. 



Figure 2. - Typical variation of amplitude ratio versus coolant wall temperature T cw showing 
singularity at T cw = T GW . 
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These programs assume that the hot-gas-side temperature is measured. Therefore, 
the coolant wall temperature is calculated from the following equation: 

If the wall temperature is measured at a point other than the hot-gas surface 
x/L ~ 0, the programs must be modified to account for the new location. 



INPUT CARDS 
Program BL 


The BL program requires the following cards for input information: 

Card 1: Three values with 312 format (table I, line 1) containing 
IRDG reading number 

ISTA station number 

NR number of readings to be processed 

Card 2: Nine values with 10E8. 4 format (table I, line 2) containing 
K conductivity 

RHO density of wall material 

CC specific heat of wall material 

L thickness of wall material 

XL ratio of distance measured from hot-gas side of wall and temperature meas- 
uring location to wall thickness 
TG mean temperature of hot gas 

TWO wall temperature, a function of time, at x = 0 

TWL wall temperature, a function of time, at x=L 

TC mean temperature of coolant 

Card 3: Contains pairs of values (table I, line 3) under 10E8. 4 format showing 
WRAD(l) frequency of temperature oscillation 

OMDTG(l) ratio of amplitude of wall temperature oscillation to amplitude of sinu- 
soidally driven hot-gas temperature 
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WRAD(NR) 

OMDTG(NH) 

The program read in HR pairs of values, up to five pairs in card 3. If more than 
five pairs are desired, more cards may be added. The NR pairs will continue to be read 
under the same 10E8. 4 format. 


TABLE I. - SAMPLE INPUT FOR PROGRAM BL 


STATEMENT | 
NUMBER 

o 

FORTRAN STATEMENT 

IDENTIFICATION 

1 2 3 U 5 

6 

7 8 

9 10 II 12 

13 16 IS 16 17 IS^ 

1 19 20 21 22 23 26 

25 26 27 28 29 Jo| 

31 32 33 36 35 36 

37 58 39 60 tl 62 

63 66 65 66 67 tel 

69 50 51 52 53 56 

55 56 57 58 59 60 

61 62 63 66 65 66 

67 6e 69 70 71 72 

75 76 75 76 77 78 79 8 

1 1 

1 














.20 7 

0 

-3 

. 2 8 

6 0 + 0 . 

1170+0 

. 6 00 0 

- 1 .16 

6 7+0 . 

9 600 + 3 

.65 80 

+3 .65 

8 0+3 . 

5 706 + 3 


.19 7 

0 

+ 0 

.15 

38 + 0 












Program JCL 

The JCL program requires the following cards for input information: 

Card 1: Three values with 312 format (table II, line 1) containing 
IRDG reading number 

ISTA station number 

NR number of readings to be processed 

Card 2: Eight values with 10E8. 4 format (table II, line 2) containing 
K conductivity 

RHO density of wall material 

CC specific heat of wall material 

L thickness of wall material 

XL ratio of distance measured from hot-gas side of wall and temperature meas- 

uring location to wall thickness 
TG mean temperature of hot gas 

TGW wall temperature, a function of distance 
TC mean temperature of the coolant 

Card 3: Contains pairs of values (table n, line 3) under 10E8. 4 format showing 
WRADA(l) frequency of temperature oscillation 

PHIRA(l) phase lag angle between driving gas temperature and responding wall 
temperature 
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WRADA(NR) 

PHIRA(NR) 

The program reads in NR pairs of values up to five pairs in card 3. If more than 
five pairs are desired, more cards may be added. The NR pairs will continue to be read 
under the same 10E8. 4 format. 


TABLE II.. - SAMPLE' INPUT FOR PROGRAM JCL 


STAItMFN r 
NUMRER 

1 

FORTRAN STATEMENT 

IDENTIFICATION 

1 2 3 b 5 

6 

7 8 9 10 ll‘ 1? ' 

13 lb IS 16 17 IS 

■ 19 20 2\ 22 25 2b 

25 26 27 28 79 3c| 

31 32 33 3b 35 56 

57 38 39 bO bl b2 

63 bb bS b6 b7 bB 


55 56 57 58 59 ooj 

61 62 63 6b 65 66 

67 68 69 70 71 72 

73 7b 75 76 77 70 79 8 

1 1 

1 








b9 50 51 52 53 5b j 





. 207 

0 

-3 .28 

60+0 . 

1170+0 

. 60 00 

-1 .16 

6 7+0 . 

96 0 0 + 3| 

. 65 8 0 ! 

+ 3 . 57 

06 + 3 



.1 9 7 

0 

+04. 3 9 

65+1 












OUTPUT 

The outputs of both BL and JCL programs display the descriptive title followed by a 
list of the input values with their program labels. The results of the calculations appear 
under the heading OUTPUT shown in the next sections. 


Program BL 

TCWA calculated coolant wall temperature 

TEST1 heat flux to coolant 

TEST2 heat flux from hot gas 

TEST ratio of heat flux to coolant to heat flux from hot gas (ideally should be 1.0) 

HCT time averaged convective heat-transfer coefficient on coolant side of wall 

HGT time averaged convective heat-transfer coefficient on hot-gas side of wall 

ERROR (1.0- TEST) which is used by the program to test against the convergence 

/ 

criterion 

R ratio of convective heat-transfer coefficients 

THETA final calculated value of ratio of the amplitude of wall temperature to amplitude 
of sinusoidally driven hot-gas temperature 
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Program JCL 


TCW calculated coolant wall temperature 

TEST1 heat flux to coolant 

TEST2 heat flux from hot gas 

TEST ratio of heat flux to coolant to heat flux from hot gas (ideally should be 1.0) 

HCT time averaged convective heat-transfer coefficient on coolant side of wall 

HGT time averaged convective heat -transfer coefficient on hot gas side of wall 

ERROR (1.0- TEST) which is used by the program to test against convergence cri- 
terion 

R ratio of convective heat-transfer coefficients 

PSIC convective heat-transfer parameter on coolant side of wall 

PSIG convective heat-transfer parameter on hot-gas side of wall 

GR check on the zero of eq. (B51) which should be equal to 0. 0000 


COMPUTER PROGRAMS 

These programs are written in FORTRAN IV and are operational on the IBM 7094-2/ 
7044 direct- coupled system of the Lewis Research Center. Program BL uses approxi- 
mately 5187 )jq storages and Program JCL uses approximately 4261 )^q storages. Flow 
charts of the programs are shown in figures 3 and 4. 

Machine running time for both programs is minimal, each using less than 1 minute 
per case. 
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PROGRAM BL - AMPLITUDE RATIO GIVEN 


/ l ✓ 

I 

2.12 1 

I 1 


/ / 

/ -READ BCD / 

/ UNIT 5 / 

/ FORMAT 100 / 

/ / 


I 2 

********************* 

* I/O LIST 16 * 

* IRDG, ISTA, NR * 
********************* 


I 3 

/ / 

/ -READ BCD / 

/ UNI T 5 / 

/ FORMAT 101 / 

/ / 


I A 

********************* 

* I/O LIST IS * 

* K, ft HQ, CCt L, * 

* XL, TB, TWO, TWL, * 

* TC * 

********************* 


I / 


/ I 


/ / 

/ READ BCD / 

/ UNIT 5 / 

/ FORMAT 101 / 

/ / 


I 6 

********************* 

* I/O LIST IS * 

* (HRADf I ) tOMDTGI l ) »* 

* 1=1, NR) * 

********************* 


I 7 

********************* 

* PRINT TITLE AND * 

* INPUT * 

********************* 



/ / 

/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 102 / 

/ / 


I 9 


/ / 

/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 103 / 

/ / 
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/ I 

I 

I 

I 18 

********************* 

* I/O LIST IS * 

* K, RH3 , CC, L, * 

* XL, TG, TWO, TWL, * 

* TC, WRAOIII, * 

* OMDTG I I ) * 

********************* 

I 

I 

I 19 

********************* 

* INITIALIZE * 

* STARTING VALUES * 

* FOR ITERATION * 

********************* 


I 

I 20 


TGW = TWO I 

I 

TCW = TC ♦ 10. I 
I 

153 = l I 


I 

1 21 
********************* 

* DEFINE HC * 

* BOUNDARY * 

* CONDITIONS AND * 

* CALCULATE R RATIO * 
********************* 


2.04 1 

9 I 22 


1 1 * I 
I ( TG-TGW J / I TCW-TC ) I 
I I 
I HG = - K * I 
I ( TCW-TGW I / ( L* ( TG- I 
I TGW) I I 
I I 
I HC = R * HG I 


* I/O LIST IS * 

* IRDG, ISTA * 
********************* 

I 

I 

I 11 


I ALPHA * I 

I K/(RHO*CCJ I 


I 12 


I BEGIN DO LOOP I 

I 2 I - 1, NR I 
I I 

l SET I * l I 



/ / 
/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 104 / 

/ / 


I 

I 14 

********************* 

* I/O LIST IS * 

* I * 

********************* 

I 

I 

1 15 


/ / 

/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 105 / 

/ / 


I 

I 16 


/ / 

/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 106 / 

/ / 


17 


I 

I 23 

********************* 

* CALCULATE AND * 

* COMPARE AMPLITUDE + 

* RATIOS * 

********************* 

I 

I 

I 24 


ETA = 

SQRT IWRADI I ) / (2.0* 
ALPHA) ) 

ETAL * ETA * L 

ANG = 2.0 * ETAL 

SANG = SIN! ANG) 

CANG = COS ( ANG ) 

EXPANG = EXP (ANG) 

THE TA ( 2 ) = 
OMDTG( I ) 

ANGl = ETAL * XL 

EXP1 = EXP1ANG1) 

ANG2 = ETAL * 
(2.0-XL) 

EXPZ = EXP ( ANG2 ) 

C0S1 = COS (ANGl) 

SIN1 = S INI ANGl ) 

PHIC = K * ETA/HC 

PHIG = R * PHIC 

PGPPC = PHIG * 
PHIC 

PGPC2 = PHIG * 
PHIC * 2.0 

F = PHIC 


I / 


/ / / 

/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 107 / 

/ / 


Figure 3. - Flow chart for program BL 



PROGRAM BL 


AMPLITUDE RATIO GIVEN 


/ 


/ I 
I 
I 

! 1 

E * 1.0 - F I 

PCP1 «= PHIC ♦ 1.0 I 

D = F * CANG ♦ I 

PCP1 * SANG I 

C = PGP1 * CANG - I 
PHIC * SANG 1 

B * (PGPPC-PGPC2 I I 
+ EXPANG * I 

( ( PGPPC+PGPC2I* I 
CANG+I PGPPC+l .0) * I 
SANG l ! 

A = UO - PGPPC * I 
EXPANG * I 

t ( PGPPC+ PGPC2 ) * ! 

SANG-4 1.0+PGPPC1 * 1 
CANG ) I 

DELI « C * A - D I 

* B 1 

GEL2 » D * A ♦ C I 

* B I 

CEL3 * F * A - E I 

* B I 

CELA = E * A ♦ F I 

* 8 I 

NUH 1 * I 

<DEU*SINI-DEL2* I 
COSH * EXP2 I 


1 

I 

I 2 


NUH 2 = 

1DEL9*C0S1-DELA* 
SIN1) * EXP1 

NUH = NUH 2 - NUMl 

DEN * 

IDEL'A*C0S1+DEL3* 
SINl) * EXP1 - 
(DEL2 *C0Sl+DEL2* 

simi # exp 2 

A2P62 = A**2 + 
B**2 

CON = 

SCRT(NUM**2+DEN** 

2) 

AMPOTG = 
CBN/A2PB2 

DIFF =« AMPDTG - 
OMOTGU) 


I 

I 3 


I I BALAN 
I 3 I 

I . I ( TC W, l GO, D4 FF ) 
I 0 I 
I 1 I 


I 

I 

I 

* A 

* * 


TRUE * 

* IGB 

I * 

I * 

I * 

I **** i 

-* 1 ♦ 

* 22 ♦ 
♦*** 

9 


I 

I 

IJK = 1 

TCHAdJKJ * TCW 
HCT * HC 
HGT = HG 

THETA (A) = AHPDTG 

TEST1 = HCT * 
(TCW- TC) 

T6ST2 = HGT * 

( TG-TGW) 

TEST = 

TEST 1/TEST2 


FORMAT STATEMENTS 


... 13 
EXIT. 


I 

I 6 

********************* 
* PRINT OUTPUT * 

********************* 

I 

I 

I 7 


/ 


/ 

/ WRITE 

BCD 

/ 

/ UNIT 

6 

/ 

/ FORMAT 

109 

/ 

/ 


/ 

I 

I 8 

/ 


/ 

/ WRITE 

BCD 

/ 

/ UNIT 

6 

/ 

/ FORMAT 

110 

/ 

/ 


/ 

I 

I 


9 

/ 


/ 

/ WRITE 

BCD 

/ 

/ UNIT 

6 

/ 

/ FORMAT 

111 

/ 

/ 


/ 


I 

I 

I 10 

********************* 

* I/O LIST IS * 

* TCHAD JK), TEST1, * 

* TEST2, TEST, HCT, * 

* HGT, ERROR, R, * 

* THETA ( A) * 

********************* 

I 

1 


I 


11 


LO/EQ* * 

* ADD 1 TO I , * 

I * COMPARE TO NR * 

I * * 

I * * 

I **** * * 

-* 1 * * 

♦ 13 * I HIGH 

**** I 


I 


12 


********************* 

* END OF DO LOOP * 

* VARIABLE I * 

********************* 

I 



. 1 . 01 . 


Figure 3. - Continued. 



PROGRAM BL - AMPLITUDE RATIO GIVEN 


SUBROUTINE BALAN 


PARAMETERS CX,K,Y) 


/ BALAN / 


X 

2.03 M 

I 1 

********************* 
* SUBROUTINE BALAN * 
********************* 
I 


I 2 

***************** 

•COMPUTEO GO TO * 

* FOR K ♦ 

* 1 AT 3.03* 

* 2 AT 3.05* 

***************** 

! 

i 


BRANCH ACCORDING TO THE 
VALUE OF VARIABLE K 


/ 1 / 

I 

3*02 “I 

I 3 

XI - X I 

Y1 - Y I 

X - X ♦ X/10. I 

K » 2 I 

I 

I 

1 

... A 

I exit! 


/ z / 


I 

3.02 1 

I 5 


I X2 = X I 

I 1 

I Y2 * Y I 


I 


* 6 
* * 

* * 


V 


ABS(Y) -LT. 
l.E-4 



I TRUE 

I 

I 



I 

I 

I 

7 

I K * 15 I 

1 

1 

I 8 

I 

X » 

I 

I 

IX2-X1 J/IY2-Y1 I 

* I 

I 

(-Y1) ♦ XI 

I 

I 

XI * X2 

I 

I 

Y1 = Y2 

I 


I 

I 

I 

... 9 

! exitI 


Figure 3. -Concluded. 


Program BL 


$ IE FTC BL 

C 

C 

C** 4 4 4 ***************************** ********** *************************** 

C CONVECTIVE HEAT-TRANSFER COEFFICIENTS FROM HARMONIC TEMPERATURE 

C OSCILLATIONS - AMPLITUDE RATIC GIVEN 

C 4 ****** *** 44 *************************************************** ******** 

c 

c 

C IRDC-=RE AD ING NUMBER 

C ISTA=STATIQN NUMBER 

C NR =NUMBER OF READINGS TO BE COMPUTED 

C (^THERMAL CONDUCTIVITY OF lnALL MATERIAL 

C RHO = CEN S ITY OF kALL MATERIAL 

C CC=SP ECIFIC HEAT OF LL MATERIAL 
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nnnfinn^nnonoooririonnononfioonn 


L = I P 1CKNE S S OF WALL MATERIAL 

XL -RA TIG OF DISTANCE MEASLREO FROM HOT GAS SIDE OF THE WALL 
TEMPERATURE MEASURING LCCATICN TO WALL THICKNESS 
TC-=NEAN TEMPERATURE CF HCT GAS 
T WO = WA UL TEMPERATURE, FLNCTICN CF TIME, X=0.0 
TWl=WALL TEMPERATURE, FUNCTICN CF TIME, X= L 
T C -ME AN TEMPERATURE OF THE CCCLANT 
wRAD =FR EQUENC Y OF TEMPERATURE CSCILLATICNS 
OMCTC-=RA1IO OF THE AMPLITUDE CF THE WALL TEMPERATURE TC THE 

AMPLITUDE OF THE SINUSOIDALLY DRIVEN HCT GAS TEMPERATURE 
AL PF A - TP ERMAL OIFFUSIVITY 

ET A = FR EQUENC Y ANO WALL PROPERTY PARAMETER 
PCA I 1 ) =LU WER LIMIT OF COCLANT SIDE CCNVECIIVE HEAT-TRANSFER 
COEFFICIENT 

PC A ( 3 ) =UP P ER LIMIT OF CCCLANT SIDE CCNVECTIVE FEAT-TRANSFER 
COEFFICIENT 

FCT = TI ME AVERAGED CONVECTIVE HEAT-TRANSFER CCEFF IC I ENT CN 
COOLANT SIDE OF THE WALL 

FG T =T I M E AVERAGED CONVECTIVE HEAT-TRANSFER CCEFFICI ENT CN HOT 
GAS SIDE OF THE WALL 
IJK ^NUMBER OF ITERATIONS 

R=RA1IG OF CONVECTIVE HE A 1- TR AN SFER CCEFFICIENTS 
TCWA=CALCULATED COOLANT WALL TEMPERATURE 
TEST 1 =FE A 1 FLLX TO THE CCCLANT 
T E ST2 = P EA T FLLX FROM THE HCT GAS 


***> 1 ******** ******** *************************************************** 
DIMENS ION TCWAC 3C0) ,WRAD(300) ,CMOTGI300) ,HCA (4) ,THETA(4) 
EQUIVALENCE (PSIC ,PHIC) ,( PSIG,PHIG) 

REAL K,L, NLM1,NUM2,NUM,NEWERR 
C * * * * ********** ********************************************************* 
C READ INPUT 

C < ********************************************************************** 
1 R EAC( 5 » 1C C ) IRDG , I STA ,NR 

REAC( 5,10 1) K , RHQ ,CC , L ,XL ,TG , TWC,TWL,TC 
REAC<S,101) ( WRAUII ) ,OMDTG(I ) »I = l, NR) 

C * X* * ***********************************************************.,******* 
C PRINT T ITLE AND INPUT 

C ************************** **************** **** ***** ******************** 
WRITE (6,102) 

WRITE(€,IC3) IRDG , I STA 
ALPFA=K/( RHO*CC ) 

DO 2 1=1, NR 
WR ITE( 6 ,ICA) I 
WRITE( 6 , ICS ) 

WRITE! ( » 1 C 6 ) 

WRITE! 6, 107) K ,RHO,CC , L , XL , TG , T VC , T WL , TC , WRAD ( I > , CMCT G< I) 

C**** ********** ***********************************$********************* 
C INITIALIZE STARTING VALLES FOR ITERATION 

C*** ******* ********************************* * 4 * ************* ***^^ ******* 
T 0 W = T W 0 
TCw=TC * 10 . 

I GO = I 

C* ********** ************************************************* 4 ** *****4,^,* 
C DEFINE PC BOUNDARY CONDITIONS AND CALCULATE R RATIO 

C * ************* ********************************************************* 
*5 R = < TC-TGW )/ ( TC W— TC ) 

PG=-K*( TCW— TGW)/(L*( TG— TG W ) ) 

F C = R * P C- 
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C* ** ********************************** ***********************4********** 
C CALCULATE AND COMPARE AMPLITLOE RATIOS 

C=M*>f *********************** ****************** ************************** 
ETA=SQPT( WRADil ) / < 2. C*ALPHA ») 

ETAL=E IA* L 
ANG=2 . (*E TAL 
SAN G = S INI ANG) 

C ANG = CC SI ANG ) 

EXPANG=EXP (ANG) 

TFETA! 2 ] = CMDTG ( 1 1 
ANG 1=ETAL*XL 
EXP 1= £ >P ( ANG 1 ) 

ANG2=E1AL*(2.C— XL ) 

EXP 2=E XP ( ANG2 ) 

COS 1=CCS( ANG 1 ) 

S IN 1= S IN ( ANG 1 ) 

PHC = K*ETA/HC 

P F I C = R <P F IC 

PGPPC=PHG*PH1C 

PGPC2=PFIG*PHIC*2.C 

F=PFIC 

E=1 .C-F 

PCP ] = P F IC -t 1.0 

C = F*CANG + PCP 1* SANG 

C=PCP 1*CANC~P PIC* SANG 

B = ( PGPPC-PGPC2 )+EXPANG*( ( PGPPC + PGPC2 1 *C AN G+ ( PGP PC+ 1 .0 ) *SANC ) 

A=1 .C-PC-PPC+EXPANG*( (PGPPC + PGPC2 )*SANG-(1 .O+PGPPC) *CANG ) 

DEL 1 = C*A-C*B 
C EL 2 = C *A *C*B 
CEL 3=f *A-E*9 
CEL A=E*A*F*B 

NUM 1= ( CEL I*SIN1-DEL2*CCSI )*EXP2 
NUM2=<CEL2*C0S1-DEL4*SINU*EXF1 
NLM=NLM2-NLM1 

DEN = ( CEL4*COS1 + OEL3*SINI)*EXP1- (OE L 1*C C SI + DEL2 *S I N 1 ) * EX P2 

A2PE2=A**2+B**2 

CON = SQRT( NLM**2*DEN**2) 

AMPCTG=CON/A2PB2 
DIFF=APPQ1G-0MDTG (I I 
CALL 6ALAN(TCl»tIGG»DIFF) 

I F ( IGC .L T .10 ) GO TO 9 
I JK = 1 

TCWAt I JK ) = TC W 
F C T = PC 
PGT=PG 

TFETAI <■ ) = AMPDTG 
TEST 1 = FCT* ( TC F- TC ) 

TEST2=FCT*I TG-TGW) 

TEST=TEST 1/TEST2 
ERROR = 3 .-TEST 

(;**** ******************************************* ************************ 
C PRINT OUTPLT 

C *********************************************************** ************ 

WR I TE ( 6, 1C9 > 

WRITE! <,11C) 

WRITE! 6, 111) TCWA( I JK ) , TE ST1 , TE S 12 ,T€ST,HCT ,HGT , ERRCR , R »TH ET A( A ) 

2 CONT INLE 
GO TO I 

C** ************ ******************** ****************** ******************* 
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r> o 


FORMA T STATEMENTS 

4 4 2 * 4 *# 4 ** 4444 444444 * 44444 44 * 4444 * 44444 ** 444 * 44 444444 ** 444444*444 4444*4 
IOC FORMAT ( 21 2 ) 

1C1 FORMAT ( 10E 8.4 ) 

1C 2 FORMAT { 12 CF ICON VEC T I VE HE A T- T P ANSFER COEFFICIENTS FROM FARMCNIC TE 
1MPERATURE OSCILLATIONS - AMPLITUDE RATIO GIVEN I 

1C 2 FORMAT ( 18FCREAOING NUMBER = ,I2,10X,18H STATION NUMEER = ,12) 

1C4 FORMAT (4HCRUN , 12) 

1 C 5 FORMAT(CHCINPLT) 

ICC FORMAT < 1HC»4X»1HK»8X,3HRHC , S X ,2HCC » 8 X ,1HL ,9X , 2HXU , 8 X, 2FTG, 7X,3HTWO 
1, 7X,3FTV»U ,7X,2HTC,7X,4H*RAD ,6>,5HCMDTG) 

IC7 FORMAT ( 1HC ,F 8 . 6 , 4F 10 . 4 , 4F 1C. 1 .2F10.4) 

ICC FORMA T(55FC ITERATIONS EXCEED 5C ) 

ICS FORMAKIHCOITPI T) 

11C FORMA T( 1HC»3X»4HTCIaA,5X»5HTEST1,5X,5HTEST2 ,6X ,4HTEST,6X,3HFCT, 7X, 3 
1HFGT , CX, 5FERR0R,7X,1HR ,7X,5H THETA) 

111 FORMAT ( 1HC,F8. ], 2F10.4,2F1C.6,3F10.4) 

112 FORMAT C55FCSOLLTION IS DIVERGING ) 

END 


$ICFTC eALAN 

SUBR0L1INE BALAN{X,K,Y) 

GO TO ( 1 , 2 ) , K 

1 X1=X 

Y 1 = Y 

X=XTX/ 10. 

K = 2 
RETURN 

2 X2 = X 

Y 2 = Y 

IF( ABS (Y ) .LT. l.E-4) K=15 
X=( X2-X1 ) /<Y2-Y1)*(-Y1)+X1 
X 1 = X2 

Y 1 = Y2 
RETURN 
END 


Program BL Output 

CONVECTIVE HAT- TRANSFER COEFFICIENTS FRCM HARMONIC TEMPERATURE OSCILLATIONS - AMPLITUDE RATIO GIVEN 


RECCING NUMBER = 1 


STATI CN 

NUMBER = 1 






RUN 1 










INPUT 










K 

RHO 

CC 

L 

> L 

TG 

TWO 

TWL TC 

NRAO 

OMDTG 

0.C0C2C7 

(•2 6 60 

0.1170 

C. 0600 

0.1667 

960.0 

658.0 

658.0 57C.6 

1.2390 

C.1538 

OUTPUT 










T CVnA 

TEST1 

TES72 

TEST 

HCT 

HGT 

ERRCR 

R T FETA 



615.0 

C.1344 

0.1344 

0.9998 

0.002775 

0.000445 

0.0002 

6 .2250 C • 1538 



*01* UNITC5 

* EOF. 






REC- CCOOC 

FlL= 00G02 
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PROGRAM JCL - PHASE LAG ANGLE GIVEN 


’ 3. 1415926 1 


**9t ***************** 

* READ INPUT * 

********************* 


READ BCD 
UNIT 5 
FORMAT 100 


I 4 

9 * 9 * 9 99 ** 9*44 ******** 

* 1/0 LIST IE * 

9 IR06, ISTA, NR * 

4 * 49 4 * ****** 9 * ******* 


I 


READ BCD 
UNIT 5 
FORMAT 101 


I 6 

* I/O LIST IS * 

* K, RHO, CC, L» * 

* XL, TG, TGW, TC * 
« 99 « 999999*********** 

I 

I 


READ BCD 
UNIT 5 
FORMAT 101 


**** ***************** 

* I/O LIST IS * 

* (WRADAI I ) tPHIRAI 1 1* 

* , 1*1, NR) * 

********************* 

I 

I 

I 9 

********************* 

* PRINT TITLE AND * 

* INPUT * 

********************* 


WRITE BCD 
UNIT 6 
FORMAT 102 


I 11 

/ / 
/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 103 / 


I 12 

**** ***•*****•••**••* 

* I/O LIST IS * 

* IROG, ISTA * 
99*4 ***************** 

1 

I 


/ WRITE BCD / 
/ UNIT 6 / 

/ FORMAT 104 / 


*** 9 * 4 *************** 


I 

I 19 

********************* 

* I/O LIST IS * 

* K, RHO, CC, L, * 

* XL, TG, TGW, TC, * 

* WRADA(JJ) , * 

* PHIRAt J J I * 

********************* 

I 

I 

I 20 

********************* 

* INITIALIZE * 

* STARTING VALUES * 

* FOR ITERATION * 

********************* 



IJK = 1 



TCWA(IJK) = TGW 



TCW = TGW 



WRAD = WRADA I J J ) 



PHI = PHIRA(JJ) 


I 

I 

I 22 

********************* 
* CALCULATE R RATIO * 

99 ******* *9 *99 **9*9*9 

3 

2 

.17 1 

I 23 


R * 

tTG-TGW) /{ TCW-TC ) 



ALPHA = 
K/(RHQ*CCJ 



X = L * XL 



ETA = 

SORT IWRAD/I2.0* 
ALPHA) ) 



ETAL * ETA * L 



ETAL2 = ETAL * 
2.0 



ETAX * ETAL * XL 



ANG1 = PHI * ETAX 



S INI = SINUNGl) 



C0S1 = COS(ANGl) 



ANG2 = ANGl - 
ETAL2 



SIN2 = SINIANG2) 



C0S2 = C0SIANG2) 



ANG3 = PHI - ETAX I 


SIN3 = SINIANG3) 

I 


C0S3 = COS ( ANG3 

I 


ANG4 * ANG3 ♦ 
ETAL2 

I 

I 


SIN4 * $1N( ANG4 ) I 


WRITE BCD 
UNIT 6 
FORMAT 105 


I 

I 17 


/ / 
/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 106 / 

/ / 
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/ / 
/ WRITE BCD / 

/ UNIT 6 / 

/ FORMAT 107 / 

/ / 


Figure 4 - Flow chart for program JCL 



PROGRAM JCt - PHASE LAG ANGLE GIVEN 


C0S4 * COS ( ANG4 ) 
EXPl *= EXPIETAL2) 


EXP 2 « 
E*PtETAL2*I2.0~XL)I 


EXP 3 * 
EXP(E0AL2*I l.O-Xl) I 
) 

A = - SIN1 ♦ 

EXP1 * COS2 + 

BXP2 * SJN3 - 

EX03 * C0S4 

b « cbsi - expi > 
SIN2 - EXP2 * 

COS-3 ♦ EXP3 * 

SIN4 

C = SCN1 ♦ EXPl - 

SINS ♦ EXP2 * 

SIN3 + EXP3 * 

SINA 

D = - COSI ♦ 

EXPa * COS2 - 

EXP2 * C0S3 + 

EXP3 * COSA 

E =■ STNl - EXPl ' 

3IN2 ♦ EXP2 * 

SINS - EXP3 * 

SINA 

I 

I 

I 

I F =* COSI + EXPl > 

I SIN2 - EXP2 * 

I COS3 - EXP3 * 

I SIN4 


1 


G * - S INI - 

expa * cos 2 * 

EXPl * SIN3 + 
EXP3 * COSA 

I 


LO/EQ* 

* ADD 1 TO KK, 

I * COMPARE TO 3 
l * * 

I * * 

I **** * * 


2 * 
06 * 


I HIGH 


* END OF 00 LOOP * 

* VARIABLE KK * 

********************* 


I RTPOLY 
I(N,AA,5Q»L.E- 
I7*UUtVVj CONVt 
IHH,8B,CCC,DD, 
I EE I 


I 11 

********************* 

* CHECK FOR * 

* POSITIVE REAL * 

* ROOT * 

********************* 


WRITE BCD 
UNIT 6 
FORMAT 112 


********************* 


********************* 

* EOLVE FOR ROOTS * 

* OF CUBIC EQUATION * 

********************* 


SET KT = 1 


4131 - 
1C*R*IC+D ) 1/ (R+I F*I 
G) ) 

Z(2) = 

1 2» 0*A+R* ( A+BJ )./( 
2.0*R*<F+GJ > 


ZIA) =. 1.0 
AAI2) * Z (3 1 
AA.I3) » ZU) 
AAiAI * ZU), 
AAU) ■ 1.000 


I 

I 

I 5 

I BEGIN 00 

LOOP I 

I 6 KK * 

1, 3 I 

I SET KK 

a 1 I 

I 

2.C7 1 

6 I 

6 

I CONV-I KK) 

=« CRTRN I 


LO/EQ* 

* ADO l TO KT, 

* COMPARE TO 3 


********************* 

* END OF DO LOOP * 

* VARIABLE KT * 

********************* 


Figure 4. - Continued. 



PROGRAM JCL - PHASE LAG ANGLE GIVEN 


KOUT -EQ. 1 


REAL2 * UUIKK ) 
REAL! = R6AL2 * R 


( IREAL2+ZI3) )* 
REAL2+Z (2)1 * 
REAL2 ♦ Zli) 


PSIC * REAL2 
PSIG * REAL 1 


WRITE BCD 
UNIT 6 
FORMAT 108 


I 13 

********************* 

* CALCULATE AND * 

* COMPARE COOLANT * 

* WALL TEMPERATURES * 
********************* 


I KOUT » KOUT * 1 1 


ADD 1 TO KT, 
' CCMfiARE TO 3 


********************* 

* END OF DO LOOP * 

* VARIABLE KT * 


KOUT .EQ. 1 


TCWAIIJK) * TGW - I 

HGT * L * I 

ITG-TGWI/K I 

I 

TCW * TCWAIIJK) I 

I 

TEST! * HCT * I 

ITCW-TC) I 

I 

TEST2 - HGT * I 

ITG-TGWI I 


ERROR - 1.0 - 
TEST 


WRITE BCD 
UNIT 6 
FORMAT 113 


I TCWAIIJK) = I 
I ITCWAI I JK-U+TCWAl I 
I IJK) )/2.0 I 

I I 

I TCW = TCWAIIJK) I 


* IJK .LT. 50 


I******************** 

* 1/0 LIST IB » 


********************* 


Figure 4. - Continued. 



PROGRAM JCL - PHASE LAG ANGLE GIVEN 


/ 4 / 


I 

3.15 -I 

I 1 


I R * I 

I (TG-TGWl/ITCW-TC) I 


I 

I 2 

********************* 
* PRINT OUTPUT * 

******** ************* 
I 
I 

I 3 


/ 


/ 

/ WRITE 

BCD 

/ 

/ UNIT 

6 

/ 

✓ GORMAT 

109 

/ 

/ 


/ 

I 

1 

I 4 

/ 


/ 

/ WRITE 

SCO 

/ 

/ UNIT 

6 

/ 

/ FORMAT 

110 

/ 

/ 


/ 

I 

I 

I 5 

/ 


/ 

/ WRITE 

BCD 

/ 

/ UNIT 

6 

/ 

/ GORMAT 

111 

/ 

/ 


/ 


I 

I 

I 6 

*••*•*••******••*«*** 

* I/O LIST IS * 

* TCWA(IJK), TESTl, * 

* TEST2, TEST, HCT, * 

* HGT, ERROR, R, * 

* PSIC, PSIG, GR ♦ 
********************* 

I 

2.20--*—! 

I 

5*7 
* * 

* * 

* * 

fcO/EQ* * 

* A 00 l TO JJ, * 

I * COMPARE TO NR * 

I * * 

I * * 

l **** * * 

-* I * * 

* 14 * I HIGH 

**** I 

1 
I 
I 

I 6 


* END OF DO LOOP * 

* VARIABLE JJ * 

********************* 

I 

I 

I 9 


FORMAT STATEMENTS 


/ 


I 

I 

I 


10 



/ / 
/ WRITE BCD / 

/ UNIT 6 / 

J GORMAT 114 / 

/ / 


I 

1 



1 


Figure 4. - Continued. 



PR0GRAH JCfc - PHASE LAG ANGLE GIVEN 


SUBROUTINE RTPOLY 


/ RTPOLY/ 



*+4*+++*+****+*+**+*+ 
* SUBROUTINE RTPOLY * 
********************* 


I 

I 

I 2 


I NSAVE = N I 

I I 

I NN * N ♦ 1 I 

I 

I 

I 3 


I BEGIN DO LOOP I 

I 1 I ° I* NN I 
I I 

I SET I * I I 


I 

A 114 


I HII) * A(I) I 


I 

I 

I 

♦ 5 


* * 

LO/EQ* * 

— * ADD 1 TO I, * 

* CONPARE TO NN * 

* * 

* * 


I HIGH 

I 

I 

I 

I 

I 6 

********************* 

* END OF DO LOOP * 

* VARIABLE I * 

********************* 

1 

I 

I 7 


I XK « F I 

I I 

I LT = 1 I 


I 8 

********************* 

* EXTRACT ZERO * 

* ROOTS * 

********************* 

I 

I 

I / 


/ 


Figure 4. - Continued. 



PROGRAM JCt - PHASE LAG ANGLE GIVEN 


SUBROUTINE RTPOLY 


» # 
(01 


u<ni • e. 
V(N9 = 0. 
CONVIN) * XK 
N * N - l 
NN * N * I 


l 

I 


* — I 

I 

I 


* * 

* 

no/-) 


07 . 
500 


-I INITIALIZE P.Q.ETC. 
I 




10 

I 

4 

I 

S » 

0. I 

I 

I 

PS « 

0. I 

I 

QS * 

0. I 

I 

PT * 

0. I 

I 

QT * 

0. I 

I 

XK =■ F I 


I 

I 

I 

« 

5 


* 

• 


• 

* 

* 


* 

(♦) * 


* (-) 

I * 


* I 

I * 


* I 

I 

* 

* I 

I **** 

* 

* V 

-* 7 * 

« 

— 

* 01 * 


(0) 

**** 



20 



12 


6 


I R * - H(2)/H(l> I 


Figure 4. - Continued. 



PR0GRAH JCt - PHASE LAG ANGLE GIVEN 


SUBROUTINE RTFDLY 


61.05 1 


NN M M 


********************* 

* ENO OF 00 LOOP * 

* VARIABLE J * 

********************* 


REPLACE X BY (l/XI 


7. 14 — * — I 


********************* 

* SCALE * 

* COEFFICIENTS OF * 

* POLYNOM'AL. ♦ 

****** «••*•******«*** 


• OABSIHI * 

* 2 )/H( 1 ) > - *{ 0 /+) 

* DABS (H(N I/HINN ) )*— 


LT - - LT 
H * NN/2 


BEGIN 00 LOOP 
36 J = 1, M 


BEGIN DO LOOP 
30 J * li NN 



LOGIlHtU)) 


ILO/EQ* 

-* ADO 1 TO J, 
* COMPARE TO M 


ILO/EQ* 

-* ADD 1 TO J, 
* CeWPARE T0 f 


********************* 

* END OF DO LOOP * 

* VARIABLE J * 

**************** ***** 


********************* 

* END OF 00 LOOP • 

* VARIABLE J * 


S * .5**(«7NN> I 


I **** * * 


BEGIN DO LOOP 
32 J * LI NN 


0 * 1 . 

P = - 2. 


1 H(J) • HIJ9 * S I 


IfcO/EQ* * 

I — « — * ADD -1 TO J, ♦ 

« COMPARE T0 NN * 


Figure 4. - Continued. 



PRBGRAM JCi - 2HASE LAG ANGLE GIVEN 


SUBROUTINE RTPOLY 


T.23— — I 


I Q o H( AN l/HIN- 11 ! 


I ilHIft 1— HIN-*21 *Q1 /HI I 


QUADRATIC EQUATION 


l 

-<5 * * 

* 03 * 

*W 


1.22 -I 


I P » HlNI/HIN-ll I 

9.11— 1 

I 

200 * 14 


1 

14 I 

U(N-l) * Q/UIN) 

VI N) - 0. 

VIN-1I - 0. 

I 

1.18 1 

:i6 i 

CONVIN) « XK 

CONVIN-11 = XK 

N * N - 2 


N 


1 


I l-> 
I 


10.24— 1 


P/Q 

' l./Q 


em * kid 

cm - Bill 


CI2) o 612) - 

cm * p 


SSS = P * p/4. - I 
Q I 

I 

SSSS - I 

DSQRT (OABS(SSS)I I 


BEGIN DO LOOP 
220 J * 1, NN 


HIJJ « BIJ) 


BJJ) * HMI * 
6HJ-1I * P • 
BIJ-2) * Q 


CIJ) * BOS! » 
CIJ-1) * P - 
CIJ-2) » Q 


UIN) * - P/2. 
UIN-1) - UIN) 
VIN) ■ SSSS 
VIN-ll - - SSSS 


* 4L0/EQI 

* ADD 1 TO J> * 1 

♦ COMPARE TO NN * 


ItC/EQ* • 

l——* ADD 1 TO J, * 
* CCKPARB Tfl NN • 


. 23 . 

216 


********************* 

* END OF DO LOOP * 

* VARIABLE J * 

********************* 


Figure 4. - Continued. 



PROGRAH JCL - PHASE LAG ANGLE GIVEN 


SUBROBTINE RTPOLY 
/ 


/ I 
I 


«**•«**»*•*****•***»* 
« END OF DO LOOP * 

* VARIABLE J * 

«*****••«•*****«•**»♦ 
I 
I 
I 


I € 1 * 

* 

I * 

I * 
I 

I «**• 
-A 9 * 
• 09 ♦ 
«*♦* 
64 


BIN) 


!!-/♦» 

t 

I 

1 

I 


HN1 * HIN) - 
BJN-2) * « 


HN1 


I !•/♦» 

I 

I 


* XK - * 

I * DAB6 IBI Nl /HN1 ) * 

I « * 

I * * 

I *«** * * 

-* 9 * * 

• 09 * I IO/-I 

**** l 

64 I 

I 
I 
I 


110 . 01 . 


* * 
!<-/♦) 


* * !♦» 

* XK - * 

♦DABS (B(N)/H(N) )* 


♦ ♦ 

* 

HO/-) 


(0/-)* * 

♦ DABSIB(N) ) ♦ 

I ♦ * 

I ♦ * 

I ♦ ♦ 

I **** * ♦ 

-♦ 10 ♦ ♦ 

♦ 01 ♦ I (+) 

♦ ♦♦♦ I 

70 I 

v I 

* — I 

64 I 9 

I B (NN) = H(NN) - I 
I BIN-1) ♦ Q I 

I 

I 

I 10 

*****♦♦♦♦♦♦♦♦♦♦*♦♦♦♦♦ 

♦ BIN+l) TEST ♦ 
♦♦*♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

I 

I 

I 

66 ♦ 11 


(♦) ♦ XK - * 

♦ DABS I B I NN ) /H I NN ) * 

I ♦ ) * 

I * ♦ 

I * * 

j **** ♦ ♦ 

-♦ 8 ♦ ♦ 

♦ 14 ♦ 

♦ ♦♦♦ 

200 


Figure 4. - Continued. 



PROGRAM JCl - PHASE LAG ANGLE GIVEN 


SUBROUTINE RTPOLY 


NEKTON 'RAPHEON ITERATION 


10.09 1 


9.05—*—l 


LINEAR EQUATION 


cm = hid 
E(1 1 => D.| 1 > 


BEGIN DC LOOP 
T2 J = 2V NN 


********************* 

* END OP DO LOOP * 

* VARIABLE J * 

****************** 1 *** 


-I 
I 

UIN) - R 

V(N1 * 0. 

CONV(N) * XK 

N * N - 1 

NN * N + 1 

1 
I 


* 07 * 
**#* 

500 


BEGIN DO LOOP 
108 J * 1, NN 


(Q + 1.) I 


10.10 1 


ss = ciN-ai * 
CIN-i) “ CJN1 
C4N-2I 


•to I ' 
l 


I **** 
-* 10 * 


p * p * 

(BIN1*CJN-U-B(NN) 
♦CIN-2) J/SS 

Q * Q * 
l-*B(N)*C(N) + 

B I A N ) *C I N-l ) ) / SS 


18 

********************* 

* END DF DO LOOP * 

* VARIABLE J * 

********************* 


Figure 4. - Concluded, 
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Program JCL 


$ I E F 1 C JCL 


44444****4 4** *************************** ****** ****4**** *********** **** 
CONVECTIVE HEAT-TRANSFER COEFFICIENTS FROM HARMONIC TEMPERATURE 
OSCILLATIONS - PHASE LAG ANGLE GIVEN 
44 4 4 * ******** ************************************************* ******** 


IROC-=REAOING NUMBER 
IS1A=STAT ION NUMBER 

NR “NUMBER OF DATA POINTS TC BE COMPUTED 
K =TH£RMAL CONDUCTIVITY OF WALL MATERIAL 
RHO =OENS ITY OF WALL MATERIAL 
CC =SP EC IF IC HEAT OF WALL MATERIAL 
L a TH ICKNESS OF WALL MATERIAL 

X L =RA T 10 OF DISTANCE MEASURED FROM HOT GAS SIDE OF THE WALL 
TEMPERATURE MEASURING LCCATICN TO WALL THICKNESS 
TG=MEAN TEMPERATURE CF HCT GAS 
TGW=WALL TEMPERATURE, FUNCTION CF OISTANCE 
TC=MEAN TEMPERATURE OF THE CCCLANT 
WRAOA=FREQUENCY OF TEMPERATURE CSCILLATICNS 

PH IRA =PHA SE LAG ANGLE BETWEEN THE DRIVING GAS TEMPERATURE AND THE 
RESPONDING WALL TEMPERATURE 
ALPFA=THERMAL DIFFUSIVITY 

ETA=FREQUENCY AND WALL PPCPERTY PARAMETER 
GR =CUB IC EQUATION CHECK CF RCCTS, SHOULD APPROXIMATE Z ERC 
FCT = TIME AVERAGED CONVECTIVE HEAT-TRANSFER CCEFF ICI ENT ON 
COOLANT SIDE OF THE WALL 

FGT *TIME AVERAGED CONVECTIVE HEAT-TRANSFER CCEFF IC I ENT CN HOT 
GAS SIDE OF THE WALL 
IJK=NLMBER OF ITERATIONS 

R =RA T 10 OF CONVECTIVE HE AT- TRANSFER COEFFICIENTS 
P S I G =CON VEC T I VE HEAT-TRANSFER PARAMETER CN THE HCT GAS SIDE OF 
THE WALL 

PSIC=CQNVECTIVE HEA T- TRANSFER PARAMETER CN THE CCCLANT SIDE CF 
TFE WALL 

TCW=CALCULATED COOLANT WALL TEMPERATURE 
TEST 1=FEAT FLLX TO THE COOLANT 
TEST2=FEAT FLLX FROM THE HOT GAS 
CONV=CON VERGENCE CRITERION USED FCR EACH ROOT 


44*4 *******4*********************************************************** 
OIMENS ION WRADA { 1 OC ) , PHI RA 1 1 OC) , TC WA 1 1 00 > ,Z (4 ) ,CCNV 14) 

COUBLE PRECISION HHI4) ,BB (4) ,CCC (4 ) ,DD (4 ) , EE { 4 ) , AA (4 > , UUC4 > , VV ( 4 ) 
REAL L »K 

DATA CPTRN/L.E-5/ 

P 1 = 3. 1* 15S26 
P 10 I8C =P I / 1 80 . 

■C 4 4 * * ****** * *********** ************ ♦* ******************** *********** * *♦* 
C REAC INPUT 

C *************************************************************** ******** 
l READ! 5 » IOC ) IRDG , I STA »NR 

REAC( 5,10 1) K,RHO,CC»L,XL,TG,TGW,TC 



READ( 5 , 1C1) ( WRAD A ( I ) , PH IRA ( I ) , I =1 , NR ) 

C * ** ******** *** ******************* 4 ************* 44** ******************** 
C PRINT TITLE AND INPUT 

C *** * ********** ************************** 4 4 ********* ***** ******* ******** 
WRITE (6, 1C2) 

WRITE(£,1C3) IRDG.ISTA 
DO 5 J J= 1 »NR 
WRITEU.1C*) JJ 
WRITE! <,1C5) 

WR ITE ( i t 1C 6 ) 

WRITE! (,1 Cl) K »RHQ»CC » L*XL *TG »TGP» TC «WPADA(JJ) » PH I P A ( J J ) 
£*********** ************************************************************ 
C INITIALIZE STARTING VALLES FCR ITERATION 

C * i(* ************************************************************ ******** 
I JK = 1 

TCWA1 I JK ) =TGW 
TCW = T G V 

WRAC= PRAGA (JJ ) 

PHI=PH IRA l JJ ) 

C *********************************************************************** 
C CALCULATE R RATIO 

C ******************************************************* ********** ****** 

2 R =! TC— TGW )/{ TCV-TC) 

ALPFA=K/(RHO*CC ) 

X=L*XL 

£ TA =SQ P T ! PR AD /(2.0* ALPHA) ) 

£TAL=E TA* L 
£ TAL 2 = £ TA L* 2 . 0 
ETAX=ETAL*XL 
ANG1 = P H + E TAX 
SIN 1*5 IN! ANG 1 ) 

COS 1=CCS! ANG 1 ) 

ANG2= ANG 1-ETAL2 
S IN 2= 5 IN ! ANG 2 ) 

CO$2 = CC>S( ANG2 ) 

ANG3*P HI- E TAX 
S IN3=S IN! ANG3 ) 

COS3=COS( ANG3 ) 

ANG4*ANG3 + ETAL 2 
S IN 4= S IN! ANG4 ) 

C0S4=CCS( ANG4) 

EXP 1= E >P ! ETAL2 ) 

EXP 2= EXP! ETAL2*! 2.0-XL) ) 

EXP 3= EXP! ETAL2*! l.fl-XDI 
A=— SIN 1 +E XP l*COS2 + E XP2* SI N3— E XP3*C0S4 
B=C0S1-EXP1*SIN2-EXP2*CCS3+EXP3*SIN4 
C=$ IN 1 *EXP 1* S IN 2+E XP2* SIN3 + EXP3*SI N4 
D=-COSHEXPl*COS2-E XP 2*C 0 S3*E XP 3*C OS4 
E = S IN 1-EXP 1*SIN2+EXP2*SIN3—EXP3*SIN4 
F=COS HEXP1*SIN2-EXP2*CCS3-EXP3*SIN4 
G=- SIN 1-E XP l*CO S2*E XP2*SI N3+E XP3*C0S4 
C * ********************************************************************** 
C SOLVE FOR ROOTS OF CUBIC EQUATION 

q *************************************************** ******************** 
Z( 3 ) = ! C*R*!C+D ))/(R*(F+G) ) 

Z( 2) = ( i.G*A+R*(A+B)l/!2.0*R*(F+G)) 

Z( 1 ) = E /( 2 .C*R*!F + G> > 

Z ( 4 1=1 .0 
AA( 21=24 3 ) 
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AAC 21 = 2(2 ) 

AA<4) = 2(1 ) 

AA( 11=1. ODC 
DO t KK = 1 » 2 

I CON V( KK ) = CRTRN 
N = 3 

CALL RTPOLY(N,AA,50,l.E-7,UUy * V , CC-NV ,HH , e B ,CCC » CD, EE ) 

Q 4444 4** ***************************** ********** *4* ************** ******** 

C CHECK FOR POSITIVE REAL ROOT 

£ | 44 44***444 44 * 4 ********************** 4 ** 44444 * 4 4 4 * 44 44 ** 4***4********** 
KOU T=C 
CO 7 K 7=1 » 3 

IF( VV( KT 1 .NE.O.DC) GO TO 7 
KK=KT 

KOL 7=KGLT + 1 
1 CON TIN IE 

IF(KOLI.NE.O) GO TO 8 
WRITE! 1 12) NR 
GO TO 5 

fi IF(KOLT.EQ.l) GO TO 11 
KOUT=C 
DO 9 K 7= 1 y 3 

IF C UU( KT 1 . LE . C . 0 C ) GO TO 9 
KK=KT 

KOUT=KCLT + 1 
9 CONT IN LE 

IF(KOLT.EQ.l) GO TO 12 
1C WRITEItyliai NR 
GO TO * 

II IF(UU( KK ) .LE.O.DC) GO TO 1C 
12 REAL2=IL( KK) 

REAL 1=REAL 2*R 

GR = ( (REAL 2*Z( 3 ) >*RE AL2 + Z ( 2 ) 1 * PE AL2 + 2 ( 1 1 
HCT = (K*ETA ) /REAL 2 
FGT=(K*ETA)/REAL1 
P S I C=R EAL 2 
P S I G=R f AL 1 

Q4 44 4 4444 *4 4 444444444*4**44*4**4*****4*4* ********* ************** ******** 

C CALCULATE AND COMPARE COOLANT WALL TEMPERATURES 

£4**4 4* 444* 44 ******** ************44************ ***** ************ ******** 

IJK=IJK+1 

TCW A I IJK)=TGW-HGT*L*(TG— TGW1/K 
TCW =TC VA ( IJK) 

TEST1=FCT* (TC W— TC ) 

TEST2=FGT*(TG-TGW) 

T£ST=TEST 1/TEST2 
ERROR= 1.0-TEST 

IF( A8 5 ( ERROR ) *L T* • Cl ) GO TC 4 
TCW A( IJK ) = (TCWA ( IJK— 1 ) + TCWA(l JK ) )/2.0 
TCW=TC lfA( IJK I 
I F < IJK .LT.5C) GO TO 2 
WRITE (CylCE) 

GO TO 5 

A R = ( TG- 1C-W J/I TCV-TC) 

C *********** *4* ****** *4*4*4**4******4**************** ******** *********** 

C PRINT OUTPUT 

Q* *4 ******** 4*4 4*4****4*************** 4*4* 4* 4***44*4 4 4*4* 44* * 44* * 44*4*** 

WR I TE ( C, 1C9) 

WR ITE ( C, 1 1C) 
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WRITE U. ill) TC WA(IJK) ,TEST1,7EST2 ,TEST ,HCT ,HGT » ERROR , R , PS IC, PSIG, 
1GR 

5 CON 1 IN l E 
WR ITEi t , 1 1A ) 

GO TO J 

CUS ******* **** 44444444444444 44444 4*444444*4*44*44 4 44 4 4 444 4 4*44444 

C FORMAT STATEMENTS 

C 4 444 4*4*4*4 444 4 4444 4 444444444 * 4444 4 4 4 44 *4 4*4 4444 44 44 4444 444 4444 44444444 
1GC FORMA T12I2) 

1C 1 FORMAT ( 10E8.A) 

102 FORMAT ( 12CH1CONVECTI VE HEAT-TRANSFER COEFFICIENTS FPCM HARMONIC TE 
IMP ERA TIRE OSCILLATIONS - PHASE LAG ANGLE GIVEN J 

1C 2 FORMAT < 18HCREAO ING NUMBER = ,I2,10X,18H STATION NUMBER = ,121 
1 C A FORMAT ( AH OR UN » 12) 

ICE FORMATUHC INPUT) 

IOC FORMAT ! 1HC»AX»1HK, 8X ,3HRHC »BX,2HCC,8X , 1HL ,9X,2HXL,8X,2FTG, 7X,3HTGW 
1, 8X,2HTC, 6X, 5H WRADA » 5X , 5HPH l R A ) 

10 1 FORMAT (1HC,F8.6,AF10.4,3F1C. 1 .2FI0.4) 

ICE FORMAT (55HCI TER ATIONS EXCEED 50 J 

ICS FORMAT (7HC0UTPUT > 

11C FORMAT! 1HC,4X,3HTC W, 5 X,5H TEST 1, 5X.5HTEST2 ,6X»4HTEST ,6X , 3HHCT , 7X , 3H 
1HGT, EX ,5H ERROR » 7X » 1HR »8X»4HPSIC ,6X ,4HPSIG ,7X,2HGR) 

111 FORMAT! 1HC,F8.1,3F10.4,2F1C.6 .5F10.41 

112 FORMAT (40HCERR0R IN SUBR RTPOLY FOR DATA POINT NO. 15) 

112 FORMAT (63HCA UNIQUE POSITIVE REAL RCCT CANNOT EE FOUND FOR DATA PO 
1INT NO . I 5 1 
1 14 FORMAT! 1H1 ) 

END 


$ I EFTC POLRT 


SUBROUTINE RTPOLYIN.A ,L,F ,U, V ,CCNV ,H ,6 ,C ,0, E) 

N ORDER OF INPUT POLYNOMIAL 

A COEFFICIENTS OF INPUT POLYNOMIAL IN CESCENCIN ORDER 

L MAXIMUM NUMBER OF ITERATION FCR EACH ROOT 

F OESIREO TOLERANCE ! REMAINDER / CONSTANT TERM) 

U REAL PART OF ROOT 

V IMAGINARY PART OF RCCT 

CON V CONVERGENCE CRITERION USED FOR EACH ROOT 

IF AFTER L ITERATIONS CONVERGENCE WITHIN F FAILS, 
THE TOLERANCE IS RELAXED BY 10. 

WORKING STORAGE 
STORAGE 
STORAGE 
STORAGE 
STORAGE 


H 
B 
C 
0 
E 

DIMENS ION A (40) ,U<40) , V(4C) ,CCNV<40) ,H(40) ,B!40) ,C(40 ) ,0(40 ), 
OOUBLE 
COUBLE 
N SAVE 


WORKING 

WORKING 

WORKING 

WORKING 


FIXED 
D .P. 
FIXED 
S.P. 
D.P. 
C.P. 
S.P. 
S.P. 
S.P. 
C.P. 
C »P « 

C. P. 

D. P. 
D.P. 

E( 40 ) 


NN = 

DO 1 
HC I ) 

XK = 

LT = 
EXTRAC 
IF! H(N 
U(N > = 
V! N ) = 
CON V( N 
N = N- 


PR EC I SION 
PR ECISICN 
= N 
*1 

= 1 »NN 
A! I ) 


A'» U » V »H , B ,C,D,E 
P,Q,R,S,PS,CS,PT,CT,HH,HN1,XK,SSS,SSSS 


1 ZERO ROOTS 
N)) A, 3, A 
C. 

C. 

) = XK 
1 


IBACCC1C 
IBACCC2C 
IB4CCC2C 
I BAOOCAO 
IBA0CC5C 
IBACOCtC 
IBACCC7C 
IB40GC8C 
IB4CCCSC 
IBACC1CC 
IBAC011C 
IBAOC 120 
IBACC130 
I B ACC 1 AO 
IBAC015C 
I84C016C 
IBAOOITC 
IB4C018C 
I B400 ISO 
IBA002CC 
IBACC21C 
IBACC22C 
IBACC23C 
IB ACC 2 AG 
IBACC25C 
IBAC0260 
IBAOC 2 70 
I84CG280 
IBACC2SC 
IBAOC 2CC 
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NN = N + 1 
GO TO ; 

4 IFI N ) fOO , 5CC, 10 
C INITIALIZE P,Q,ETC. 

I( s - c. 

PS = c. 

QS = C. 

PT = C. 

QT = 0. 

XK = F 

IF IN-1) ECO, 12,20 
12 R = -H2I/H1) 

GO TO ICO 
2C NN = N+l 

C SCALE COEFFICIENTS OF POLYNOMIAL. 

1=0 

DO 3C J = 1,NN 
I F( HI J I) 24,30, 24 
24 1= I+L0G2I HUH 

3C CONTINUE 
FN = N 1 
S=.5**l I/NNI 
DO 22 J = 1,NN 
22 H(J )=H ( J >*S 

IF( DAB S(H( 21/HI 111 - DAB S IH I N) /H II 
C REPLACE X BY I 1/X) 

34 LT = - LT 
M = NN/2 
DO 26 J = 1, M 
HH = H(J) 

MMM =N+2-J 
HIJ) = HI MMM I 
36 HI MMM ) = HH 
4 C IFIQS) l Z ,44,42 
42 P = PS 
Q = QS 
GO TO 54 

44 IFIHN-DI 46,46,48 
46 Q = 1. 

P = -2. 

GO TO EO 

46 Q = HI NN ) /HI N- 1 1 
S9S P = I H(N)-HIN-2)*Q)/HIN-1) 

5C I FI 2— N ) 52,199,4 
52 R = C. 

C BEGIN ITERATIONS 

54 II = C 

C BAIRSTCfc ITERATION 

6 C BI 1 ) = HI 1 J 

cm = el n 

3(2) = H ( 2 ) — B 1 1 >*P 
C I 2 I = BIZ) - CI1I*P 
00 62 J = 2,NN 

BIJ) = FIJI - B I J— 1 )*P - BI J-2) *C 
62 CIJ) = BIJ) - C(J-i)*P - C I J-2 ) * C 
IFt BIN 1 ) 162,64,162 
162 HN1 = MNI - BIN-2) *Q 
IFIHN1 ) 164, 166,164 
164 IFIXK - DABS! BIN1/HN1 ) ) 7C,7C,64 
166 IFI HI N ) ) 168.17C, 168 
166 IFIXK - DABS(B(N)/HIN) )) 7C,7C,64 
17 C IFI CABSIB IN ))) 7C, 70,64 
64 BINN) = HINN) - BIN-1)*Q 
C BIN 41) TEST 


IB4C0210 
IB400320 
IB4CC23C 
I B4CC340 
IB4CC35C 
IB4CC260 
IB4C027C 
IB40C28C 
IB4CC2SC 
IB40C4CC 
IB40041C 
IB40042C 
IB4CC430 
IB4CC440 
IB4CC450 
I B40C46C 
IB4C047C 
I84C048C 
IB4GC490 
IB40G50C 
IB4C0510 
IB4CC520 
IB40052C 
IB4CC54C 

)) 34,40,40 IB4CC55G 

IB4CC56C 
IB4CC57C 
1B4C0580 
IB4CC5SC 
IB4C06CC 
IB4CC61C 
IB4CC62C 
IB40C62C 
IB4CC64C 
IB4CC65C 
IB4CC66C 
I B40C67C 
IB4CC68C 
I B4CC6SC 
IB4CC7CC 
IB40C71C 
IB4CC72C 
IB4CC72C 
I840C740 
IB40G75C 
IB4C076C 
IB40C77C 
IB4CC78C 
I840C79C 
I B4CC8CC 
I B4CC81C 
IB4CC82C 
IB4CC62C 
IB4CC84C 
IB4CC65C 
I B40086C 
IB40C67C 
IB40C88C 
IB4C085C 
IB4Q09CC 
IB4C051C 
IB4CC920 
I B4C093C 
I B4CC54C 
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66 

IFIXK - DABSIB(NN)/HINN) 1 ) 70 ,70,200 

IB4CC95C 


NEWTON 'RAPHSON ITERATION 

IB40096C 

7 C 

0( 1 ) = HI 1> 

IB40C970 


Ell) = Dll) 

IB40C98C 


00 72 J = 2 S NN 

IB40C990 


C( J ) = HU) + DIJ-1)*R 

IB4C1CCC 

72 

E( J ) = Cl J ) + E (J- 1)*R 

18 40 1C 10 


IFIXK - DABS(D(NN)/H(NN)) ) 74,74,100 

IB4C1C2C 

74 

CIN) = -C (N- 1 )*P - C IN— 2)*Q 

IB4C1C3G 


SS = C IN- 1 )*CIN-1 ) - C(N)*ClN-2) 

IB4C1C40 


IFISS) 16,78,76 

IB401050 

16 

P = P + I B(N)*C(N-l)-B«NN)*C(K-2))/SS 

I B401C6C 


Q = Q * I — B{ N )*C IN ) + BINN)*C(N-1) )/SS 

IB4C1C7C 


GO TO 80 

IB401C8C 

78 

P = P - 2. 

IB4G1C9G 


q = e*(Q ♦ i.) 

IB4C11CC 

ec 

IFIEIN)) 84, 82,84 

IB4C1110 

82 

R = R - 1 . 

18401 120 


GO TO 86 

IB40113C 

64 

R * R - D (NN ) /E IN ) 

IB401140 

86 

II = II + 1 

IB4C115C 


IF! II - L ) 60,90,90 

IB4C116C 

90 

PS = PT 

IB4C117C 


QS = Q 1 

IB4C118C 


PT = P 

IB4C119C 


QT = Q 

IB4C12CC 

92 

XK = 1 C •* XK 

IB4C1210 


GO TO 24 

IB40122C 


LINEAR EOLATION 

IB4C123C 

IOC 

IFILT) 102,104,104 

IB4C124C 

] C 2 

R = l./R 

IB40125C 

1C4 

U( N ) = R 

IB40126C 


VCN) = 0. 

IB4C127C 


CONVIN ) = XK 

I B4C12BC 


N = N-l 

IB 40 12 90 


NN = N + 1 

I84C12CC 


IF(N) 500 » SCO, 1C6 

1B401210 

1C 6 

DO 108 J = 1, NN 

IB 40 13 20 

1C 8 

HI J ) = D! J) 

IB4C132C 


GO TO 10 

I B 40 1 3 4 C 


QUADRATIC EQUATION 

IB4C1350 

199 

P = HIN )/H(N- 1) 

IB4C136C 

20C 

IFILT) 202,204, 2C4 

IB4C137C 

2C2 

P = P/C 

IB4C138C 


Q = l./Q 

IB4C1390 

204 

SSS = P*P/4. - Q 

IB4C14CC 


SSSS = DSQRTIDABSISSS) ) 

IB4C141C 


IF! SSS ) 206,208,208 

IB40142C 

206 

UIN) = - P/2. 

IB4C143C 


UIN-l) = LIN) 

IB4C1440 


VIN) = SSSS 

IB4C145C 


VI N-l ) = - SSSS 

IB4C146C 


GO TO 216 

IB40147C 

2C8 

IFIP) 210,212,212 

IB4C14eC 

2 1C 

U( N ) = - P/2. + SSSS 

IB401490 


GO TO 214 

I B 40 1 5 C 0 

2 12 

UIN) = - P/2. - SSSS 

IB40151C 

21-1 

UIN-l ) * Q/UIN) 

I84C152C 


VIN ) * 0. 

IB4C153C 


VI N-l ) = 0. 

IB4C154C 

216 

CONVIN 3 = XK 

IB4C155C 


CON VI N— 1 ) = XK 

IB40156C 


N = N-2 

IB4C157C 


NN * N +1 

IB 40 1580 


1FIN) 5C0 , 500 ,218 

IB40159C 

218 

DO 220 J = 1,NN 

I B4016CC 



22C H(J I » 6 ( J ) IBAClfclC 

60 TO 10 1B401E2C 

C RESTORE N, RETURN. IBAC163C 

EOC N = N SAVE IBAC16AC 

RETURN IBAC16EC 

END I84C1EEC 


Program JCL Output 


CONVECTIVE HEAT-TRANSFER COEFFICIENTS FRCM HARMONIC TEMPERATURE OSCILLATIONS - 
READING NUMBER « 1 STATION NUMBER = 1 

RUN 1 
INPUT 


K 

RHO 

cc 

L 

XL 

TG 

T GW 

TO 

0 .CCCZ07 

C .2 860 

0.1170 

C.C600 

0.1667 

960.0 

658.0 

570.6 

OLTPLT 








TCW 

TEST1 

TEST2 

TEST 

HCT 

EGT 

ERROR 

R 

624.6 

(.1146 

0.1152 

0.9949 

0.002121 

0.000381 

0.0051 

5.59C3 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, July 13, 1970, 

129-04. 


PHASE LAG ANGLE GIVEN 


WRACA PHIRA 

1.226C C . 7673 

PS1C PSIG GR 

C.5764 5*4305 -0.0000 
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APPENDIX - SYMBOLS 


Machine 

Engineering 


code 

symbol 


ALPHA 

a 

2 2 

thermal diffusivity, K/pC, in. /sec; m /sec 

CC 

c 

specific heat of wall material, Btu/(lbm)(°R); J/(kg)(K) 

ERROR 

— 

1.0 - (TEST1)/(TEST2) 

ETA 

V 

frequency and wall property parameter, 1/in.; 1/m 

GR 

— 

cubic equation check of roots, equals zero for correct solution 

HCA(l) 

— 

lower limits of the coolant-side convective heat-transfer coef- 
ficient 

HCA(3) 

— 

upper limits of the coolant-side convective heat-transfer coef- 
ficient 

HCT 

h c 

time averaged convective heat-transfer coefficient on the coolant 
side of the wall, Btu/(in. 2 )(sec)(°R); W/(m 2 )(K) 

HGT 

h G 

time averaged convective heat-transfer coefficient on hot-gas side 
of wall, Btu/(in. 2 )(sec)(°R); W/(m 2 )(K) 

IJK 

— 

number of iterations 

IRDG 

— 

reading number 

ISTA 

— 

station number 

K 

K 

thermal conductivity of wall material, Btu/(in. )(sec)(°R); 
J/(m)(sec)(K) 

L 

L 

thickness of wall material, in. ; m 

NR 

— 

number of readings to be computed 

OMDTG 

0/at g 

ratio of amplitude of wall temperature to amplitude of sinusoidally 
driving hot-gas temperature 

PHIRA 

<P 

phase-lag angle between driving gas temperature and responding 
wall temperature, rad 

R 

R 

ratio of convective heat-transfer coefficients, h c /h G 

REAL1 

^G 

convective heat -transfer parameter on hot-gas side of wall (JCL 
program), Kr]/h G 



Machine Engineering 
code symbol 


REAL2 

if/ 
r c 

convective heat-transfer parameter on coolant side of the wall, 
Kl/h c 

RHO 

p 

3 3 

density of wall material, Ibm/in. ; g/m 

XL 

x/L 

ratio of distance measured from the hot-gas side of wall to wall 
temperature measuring location to the wall thickness 

TC 

T c 

mean temperature of the coolant, °R; K 

TCW 

T cw 

calculated coolant wall temperature (JCL program), °R; K 

TCWA 

T cw 

calculated coolant wall temperature (BL program), °R; K 

TEST1 

% 

2 2 
heat flux to coolant, Btu/(in. )(sec); W/m 

TEST2 


2 2 

heat flux from hot gas, Btu/(in. )(sec); W/m 

TG 

t g 

mean temperature of hot gas, °R; K 

TGW 

T(X, T) 

wall temperature, function of distance (measured into wall from 
hot surface) and time (JCL program), °R; K 

TWL 

T(L, r) 

wall temperature, function of time at x = L (BL program), °R; K 

TWO 

T(0, t) 

wall temperature, function of time at x = 0 (BL program), °R; K 

WRAD 

f 

frequency of temperature oscillations (BL program), rad/sec 

WRADA 

f 

frequency of temperature oscillations (JCL program), rad/sec 
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